function y = powermod(z, m, CONSTS)

    plot_data = false;
    
    E_phi_m = (func_E_phi_m(z, m, CONSTS));
    j_m_conj = conj(func_j_m_for_power_mod(z, m, CONSTS));

    y = j_m_conj.*E_phi_m;
    
    if(plot_data)
        figure; plot(z, real(y), 'b-', z, imag(y), 'r-');
        title('power_{mod} for integration');
    end

end

function y = func_j_m_for_power_mod(z, m, CONSTS)

    d = CONSTS.d;
    coef_Im = (1/pi)*func_coef_Im(m, CONSTS);

    if(m==0)
        coef_Im = (1/pi)*func_a0(CONSTS);
    end
    
    y = coef_Im./((d^2-z.^2).^0.5);

end

function y = func_E_phi_m(z, m, CONSTS)

    plot_data = false;

    k0 = CONSTS.k0;
    a = CONSTS.a;
    d = CONSTS.d;

    right_conditions = check_boundary_conditions_for_modes(m, CONSTS);

    N_mn = func_norm(m, CONSTS);
    
    if(right_conditions)
        [foo, p_mn] = func_coef_of_fields(m, CONSTS);
        [foo,   E_phi_m_border,   foo]  = func_E_inside (1, m, CONSTS);
    else
        p_mn = 0;
        E_phi_m_border = 0;
    end

    coef_Im = (1/pi)*func_coef_Im(m, CONSTS);

    if(m==0)
        coef_Im = (1/pi)*func_a0(CONSTS);
    end

    y = coef_Im*((2*pi*a)/N_mn)*((E_phi_m_border)^2)* ...
        (exp(-1i*k0*p_mn*z)).*pi.*besselj(0, k0*p_mn*d);
    
    if(plot_data)
        figure; plot(z, real(y), 'b-', z, imag(y), 'r-');
        title('E_{phi_m}');
    end
        
end



